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Coordinated Spatial Pattern Formation in 
Biomolecnlar Commnnication Networks 

Yutaka Hori, Hiroki Miyazako, Soichiro Kumagai, and Shinji Kara 


Abstract 

This paper proposes a control theoretic framework to model and an¬ 
alyze the self-organized pattern formation of molecular concentrations in 
biomolecnlar communication networks, emerging applications in synthetic 
biology. In biomolecnlar commnnication networks, bio-nanomachines, or 
biological cells, communicate with each other using a cell-to-cell commu¬ 
nication mechanism mediated by a diffusible signaling molecule, thereby 
the dynamics of molecular concentrations are approximately modeled as a 
reaction-diffusion system with a single diffuser. We first introduce a feed¬ 
back model representation of the reaction-diffusion system and provide a 
systematic local stability/instability analysis tool using the root locus of 
the feedback system. The instability analysis then allows us to analyt¬ 
ically derive the conditions for the self-organized spatial pattern forma¬ 
tion, or Turing pattern formation, of the bio-nanomachines. We propose a 
novel synthetic biocircuit motif called activator-repressor-diffuser system 
and show that it is one of the minimum biomolecular circuits that admit 
self-organized patterns over cell population. 


D 


1 Introduction 

Designing the coordinated dynamical behavior of cell populations is one of the 
important challenges in synthetic biology. The study of cell population control 
is expected to provide not only the better understanding of biology but also the 
ability to build biomimetic nanomachines, or bio-nanomachines, for bioengineer¬ 
ing applications such as tissue engineering and targeted drug delivery [^[^. In 
the last decade, researchers designed a variety of biochemical circuits that per¬ 
form population-level behavior such as synchronized oscillations , population 
control [^1^, and multicellular computation [6j|^. In these works, intercellular 
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coordination was enabled by a cell-to-cell communication mechanism called quo¬ 
rum sensing, where cells communicate with each other by releasing and receiving 
diffusible signaling molecules. 


One of the earliest synthetic biocircuits for spatial pattern formation of cell 
population was demonstrated in Basu et al. |^, where sender cells transmit a 
quorum sensing autoinducer molecule and receiver cells respond by expressing 
fluorescent protein to form predefined spatial patterns. This work was later 
followed by Liu et al. 10 , where another biochemical circuit was developed to 
form a periodic pattern of cell density in a self-organized manner. From an 
engineering point of view, self-organized pattern formation is preferable since 
it does not require extra design steps of external inputs and spatial allocation 
of cells. Moreover, the study of self-organization is important in developmental 
biology to understand morphogenesis 


11 12 


Recently, Hsia et al. 13 introduced a dynamical model of biomolecular com¬ 
munication networks using reaction-diffusion equations, where diffusion-based 


cell-to-cell communication was modeled based on the Pick’s laws of diffusion 14 


It was then shown that a biochemical communication network can form a self- 
organized concentration gradient of chemical species over cell populations de¬ 
spite the averaging effect of molecular diffusion. Such self-organized spatio- 
temporal pattern is known as Turing pattern 


15 


Although the theory of Turing pattern formation was extensively studied 
over the last 60 years |14l|16] , many existing works of Turing pattern formation, 


which are often applied to developmental biology (see 17 for example), are not 
directly applicable to today’s engineering applications in biomolecular commu¬ 
nication networks because of the difference of the number of diffusible molecules. 
Specifically, in biomolecular communication networks, only a single molecular 
species, or a small signaling molecule, can diffuse between cells and most pro¬ 
teins and mRNAs are localized in a cell, while the classical works assume that 
all molecular species can diffuse in the medium. Hence, it is desirable to develop 
a novel theoretical framework that specifically targets the systematic modeling 
and analysis of the spatio-temporal dynamics of biomolecular communication 
networks with a single diffusible molecule. 


In this paper, we consider a class of reaction-diffusion systems that approxi¬ 
mately models biomolecular communication networks, emerging applications in 
synthetic biology. Our interest here is to study the structural properties of the 
systems that come from the physical constraints of biomolecular communication 
and utilize them to derive an analysis framework rather than to study general 
reaction-diffusion systems. Specifically, the goals of this paper are (i) to present 
a stability analysis tool for biomolecular communication network systems and 
(ii) to use it to study the structures of biocircuits that are required for Tur¬ 
ing pattern formation over the population of cells. From a synthetic biology 
viewpoint, it is useful to characterize the required structures of biocircuits as it 
specifies the design space of possible combinations of genetic parts. 


We first present that the local stability analysis of the reaction-diffusion sys¬ 


tems boils down to a simple graphical test using the root locus method 18 
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We show that the root locus is characterized by the dynamics of biocircuits in¬ 
side cells and the complete circuit including diffusers, which allows us to provide 
physical interpretations of our theoretical results. Using the stability/instability 
analysis tool, we study conditions for Turing pattern formation in synthetic 
biomolecular communication networks. In particular, we mathematically prove 
that the minimum biocircuit producing coordinated spatial patterns is composed 
of three genes, and moreover, certain activation-repression network patterns are 
necessary. As an example of the three-component biomolecular circuits, we pro¬ 
pose activator-repressor-diffusor motif and show that a synthetic biocircuit with 
this motif can admit Turing instability using illustrative numerical simulations. 


The organization of this paper is as follows. In Section we provide a 
control theoretic modeling framework for biomolecular communication networks 
and introduce a useful decomposition that simplifies the stability analysis of the 
system. Then, in Sectionj^ we provide a systematic stability analysis tool based 
on the root locus method and demonstrate on a novel activator-repressor-diffusor 
biomolecular circuit. Section]^ is devoted to deriving the condition for Turing 
pattern formation. Then, in Section the class of reaction networks that can 
admit Turing patterns is characterized. Finally, Section concludes the paper. 

The authors’ conference papers [T^[^ illustrated the outlines of the Turing 
instability analysis method shown in this paper, and the proofs were partly 
presented in the paper written in Japanese 21 . In contrast with these papers, 
the present paper emphasizes emerging applications in synthetic biology. For 
this purpose, all results are presented in a two dimensional spatial domain. 
In particular, a novel activator-repressor-diffuser biomolecular circuit motif is 
proposed as an example of minimum biocircuits to achieve spatially coordinated 
pattern formation. Biological relevance and challenges for implementation are 
also discussed. In addition, we provide a formal mathematical statement of our 
key result, a necessary and sufficient exponential stability condition (Proposition 
I), for the first time with a complete proof, then the relation with the existing 
theoretical result 22 is discussed. Other mathematical results are also shown 


with an updated rigorous assumption (Assumption I) and complete proofs. 


The following notations are used throughout this paper: Nq := {0,1, 2, • • • }. 
Mo+ := {r G K I r > 0}. C+ := {c G C | Re[c] > 0}. Cq-i- := {c G C | Re[c] > 0}. 
C_ := {c G C I Re[c] < 0}. /„ is an n by n identity matrix. 


2 Modeling of Biomolecular Communication Net¬ 
works 

2.1 Model description 

We consider molecular communication networks in a two dimensional medium 
U := Si X S 2 = [0, Li] X [ 0 ,^ 2 ] where a large number of cells communicate 
using a single diffusible transmitter molecule as illustrated in Fig. Each cell 
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M-n- transmitter molecule 


Figure 1: Molecular communication network. The chemical concentrations in¬ 
side cells are approximately modeled by the continuous gradient of concentra¬ 
tions a:(4, t). 


contains a genetically identical gene regulatory circuit that serves as a functional 
module such as bistable switch [^, oscillator [^, logic gates [^, and so on. 
The circuit is composed of n—1 molecular species ' ,^n-i and is 

called internal circuit (see Fig. [^. The molecules A^ 2 , • • ■ ,-Adn-i are, for 
example, transcription factors and associated mRNAs. In order for the cells 
to communicate, cells produce a transmitter molecule and release it to the 
medium. An example of such transmitter molecule is N-Acyl homoserine lactone 
(AHL), which is a small molecule that can diffuse through cell membrane and 
is used for quorum sensing of gram-negative bacteria. It should be noted that, 
in engineering applications, one could use engineered bio-nanomachines such as 
vesicles instead of living cells. 

When the cells are densely populated in the medium, the chemical levels 
in the cells can be approximately modeled as the continuous gradient of the 
molecular concentrations inside cells. Let denote the concentration of 

Mi {i = 1,2,-•• ,n) at position | := [^ 1 ,^ 2 ]^ € fl and at time t, and define 
■= - ■ ■ ,Xn{^,t)]'^. The dynamics of the concentrations 

of Ml, M 2 , ■ ■ ■ ,Mn are then modeled by 


dx{^,t) 

dt 


f{x{^,t)) -h D 


(&^x{^,t) d'^x{^,t)\ 

1 dei ^ de2 J 


( 1 ) 


where /(•) is a function representing the rate of production and degradation 
of the molecules. The matrix D is defined as D := diag(0, • • • ,0,/r) G Kg^” 
with a diffusion coefficient It should be noted that only the (n,n)-th entry 
is non-zero, since only the transmitter molecule Mn can diffuse between cells. 
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We consider the Neumann boundary condition 


n • Va;(|, t) =0; S dil, 


(2) 


where n is a vector normal to the boundary dfl. Throughout this paper, 
we assume that the model Q, the boundary condition and initial values are 
well-posed (see 26 and Chapter 14 of 1^) and that the Jacobian of /(•) is 


well-defined at x satisfying f{x) = 0. 

Let X G Ko_|_ denote an equilibrium state of a single cell, i.e., f{x) = 0. It 
follows that the right-hand side of Q becomes zero when x{^,t) = x for all 
^ G n, because the concentration of each molecular species is uniform over the 
space and the spatial derivatives in Q become zero. This implies that a; is a 
spatially homogeneous, or spatially uniform, equilibrium state of the system 0- 
It should be noted that the spatially homogeneous equilibrium state might not 
be unique. The local dynamics around x are modeled by the following linearized 
model of Q. 


dx{^,t) 

dt 


= Ax{$„t) -I- D 




de. 


+ 


5^1 


(3) 


where A G is the Jacobian of /(•) at x, and x{^, t) is defined by x{^, t) := 
x{^,t) ~ Note that the system 0 is an infinite dimensional linear time- 
invariant system. 

In Section we provide a systematic analysis tool for studying stability of 
the equilibrium point x using the linearized model (|^ . We then show that cells 
can form spatially inhomogeneous gradient of concentrations over the population 
when the system ([^ is locally unstable and satisfies certain conditions. 


2.2 Control theoretic formulation of biomolecular commu¬ 
nication networks 

In this section, we introduce a control theoretic model representation of the 
biomolecular communication network system. We show that the infinite di¬ 
mensional system (|^ can be decomposed into an infinite number of finite di¬ 
mensional single-input single-output (SISO) subsystems that account for the 
dynamics of individual spatial modes, which will be a key result in developing 
a systematic stability analysis tool in the next section. 

The linearized model ^ can be expressed as 

= Ax{^,t) + enu{^,t) (4) 

y{i,t) = elx{i,t), ( 5 ) 

where u{^,t) ■= with the Laplace operator 


:= 


92 92 


( 6 ) 
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h{s) 



Figure 2: (Left) Feedback system representation of biomolecular communication 
network. (Right) Decomposed subsystem 


and Bn '■= [0, • • • ,0,1]^ S K". The equations Q and § imply that for each 
fixed position the dynamics of the local reactions inside a single cell can be 
modeled as a SISO system with the input u{^,t) and the output Note 

that y{^,t) is the concentration of the transmitter molecule Mn at position 
For each the transfer function from u to y is obtained as 

h{s) := e^{sln - A)“^e„ (=: n{s)/d{s)), (7) 


where we define n{s) and d{s) as the numerator and the denominator of h{s), 
respectively. As a result, the system can be illustrated as the feedback system 
shown in Fig. (Left), where I is an identity operator. Note that the feedback 
system can be viewed as a multi-agent dynamical system, where the homo¬ 
geneous agents h{s) communicate with the nearest neighbors by the Laplace 
operator The stability analysis methods of a finite dimensional version 

of such systems were established in 28,29 . Here we use the same approach to 
analyzing the stability of the infinite dimensional system 


Assumption 1. We assume (A, e„) is stabilizable, (A, e^) is detectable, and 
h{s) does not have zeros on the imaginary axis, that is, n{juj) ^ 0 /or alluj gM.. 


Stabilizability and detectability guarantees the stability of the internal states, or 
molecular concentrations, that one cannot aware of and/or influence, if any. The 
readers are referred to Chapters 14 and 16 of for the mathematical details. 
We also note that the non-existence of the zeros of h{s) on the imaginary axis 
is not restrictive in many biological applications as shown after Lemma 1 in 
Section |3l 

The stability analysis of the system ([^ is not necessarily straightforward 
due to the infinite dimensionality of the system. In what follows, we show that 
the infinite dimensional system can be decomposed into an infinite number of 
mathematically tractable finite dimensional subsystems. The Laplace operator 
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can be written as 




92 02 

M^W2 


)X-= 


( 8 ) 


with the tensor product ® of the operators and the identity operators 1=^ on 
Si = [0, Li] [i = 1,2). The eigenvalues and the associated eigenfunctions 
(t>kAi) of -AtV2 are 


Xk,i 




'■= cos 




(9) 

( 10 ) 


with fc = 0,1,2, • • • and ^ = 0,1,2, • • •. This implies that we can diagonalize the 
Laplace operator by (J'(8).L')V2(L'“i(g)L'~i) with the Fourier transform operator 
T. Since (g)L'“^)(lH 2 ® ^(s)21hi+ h(s)lH 2 J(L'(g)L') = 1^.2 ® h{s)lB.i + 
h{s)lE .2 ® , the feedback system in Fig. (Left) can be decomposed into an 

infinite number of subsystems "Sk/ {k = 0,1, 2, • • • , ^ = 0,1,2, • • •) depicted in 
Fi g. (Right). The readers are referred to, for example. Chapters 4 and 10 
of [31| for the detailed derivation of the eigenvalues and the eigenfunctions. 


This decomposition corresponds to the spatial modal decomposition of the 
molecular communication network system, which is essentially the method of 
the separation of variables (see 32 and Chapter 4 of 1^). Thus, the dynamics 


of each subsystem Y,k,i represent those of each spatial mode 4'kA^)^ ^cid it 
follows that 


WkA't)(l>kA^)^ ( 11 ) 

k=0 i=0 

where is the state of the subsystem 

This intuitively suggests that the concentrations perturbed around the ho¬ 
mogeneous equilibrium x eventually come back to S, when the growth rates 
of all the spatial modes are negative. On the other hand, if the growth rate 
of a non-zero spatial mode, fpkAA with fc -|- £ > 1, is positive, we expect the 
formation of a spatially periodic pattern of concentration gradient, or a Tur¬ 
ing pattern. Therefore, the stability analysis of the molecular communication 
network is reduced to checking the stability of the finite-order subsystems 


Remark 1. In general, the eigenvalues and the eigenfunctions defined in 


and (10) can be written as 


N 


Afc — // ^ ^ 


~T~ 


N 


AAi) = 


cos 


kiT^^i 

L, 


( 12 ) 


when the spatial domain is N dimension and Neumann boundary condition is 
imposed, where k is the set of subindices (fci, k 2 , • • • , kA concerning the spatial 
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dimension and Li is the length of the i-th dimension. The pair of eigenvalues and 
eigenfunctions depends on the boundary condition. For example, the eigenvalues 
for the Dirichlet boundary condition are given by (12) with fc + £ > 1. In 
other words, the zero eigenvalue Aq o = 0 is not an eigenvalue for the Dirichlet 
condition. The associated eigenfunctions are {kiir^i/Li). We 

refer the readers to Chapters 10 and 11 of for other boundary conditions. In 
the rest of this paper, we focus on the case of two dimensional spatial domain 
D with the Neumann boundary condition, since biomolecular communication 
networks would be implemented on a two dimensional surface. 


3 Stability analysis method for linearized model 

In this section, we provide a systematic local stability analysis tool for the lin¬ 
earized model ([^ of the molecular communication networks using a control the¬ 
oretic tool. The method is demonstrated on a novel activator-repressor-diffuser 
system. In particular, we show that the instability counterpart of the method 
is useful for the design and analysis of coordinated spatial pattern formation of 
the synthetic biocircuits. 


3.1 Graphical stability analysis method based on root lo¬ 
cus 


Suppose the internal circuit inside a cell is composed of n > 2 molecular species, 
and let the matrix A G be partitioned as 


A = 


^circ I ^cd 

^dc ' ^diff 


(13) 


with Aeirc e e and Odiff < 0. Then, 

the matrix Adrc represents the structure of the internal circuit, and adifi(< 
0) is the degradation rate of the single diffusible molecule. The vectors Ocd 
and Ode correspond to the interactions between the internal circuit molecules 
Ml, M2, ■ ■ ■ ,Mn-i and the diffuser Mn (see Fig. [^. When n = 1, we define 

A — Udiff. 

Let a polynomial p\(s) be defined by 


p}.(s) := d(s) + An(s) (14) 

with d(s) and n(s) in ([^. The characteristic polynomial of each subsystem Sk,£ 
is then written as p\^ ^(s) = d(s) + Xk,en{s). The following lemma provides a 
physical interpretation of the poles and zeros of h{s) using A and Adrc- 


Lemma 1. Consider the molecular communication network system Then, 
the characteristic polynomial p\{s) of the system can be written as 


P),{s) = d{s) + Xn{s) = \sln - ^1 + A|sJ„_i - ^cird 


(15) 






where we define |sJo — ^circl = 1 when n = 1. 

The proof can be found in Appendix. The lemma states that the zeros of 
h{s) are determined from the dynamics of the internal circuit specified by Aciro 
while the poles are determined from the overall reaction dynamics including the 
diffuser Ain- Regarding Assumption 1, we note that n{s) = |s/n_i — Acird = 0 
does not have roots on the imaginary axis for almost all Acirc when all diagonal 
entries of Adrc are non-zero. This is because the coefficient of the term is 
given by Tr(Acirc)- In fact, the diagonal entries of Adrc are non-zero in many 
biomolecular systems due to degradation and autoregulation. 

It follows that the subsystem is asymptotically stable, if and only if 
all the roots of p\^ ^(s) = 0 lie in the open left-half complex plane, also called 
Hurwitz. Since each subsystem represents the dynamics of each spatial mode 
as shown in the previous section, the stability of the molecular communication 
network ([^ is guaranteed, if and only if all of the infinite number of finite 
dimensional subsystems (fc = 0,1, 2, • • • , £ = 0,1, 2, • • •) are stable. 

Proposition 1. Suppose Assumption 1 holds. The molecular communication 
network ^ is exponentially stable if and only if px^ ^^s) is Hurwitz for all 
fc,£ G No. 


Proof. We show the proof for the case of one-dimensional spatial domain in 
order to avoid notational clutter. However, generalization to multiple spatial 
dimensions is straightforward. Let denote the roots of PXk{s) = 0. 

Then, it follows that 

n 

= '^Ck,ie'^’‘-'*WkiO), (16) 


where Ck.i (i = 1, 2, • • • , n) are some constant. 

(^) Suppose all of the polynomials Px^{s) = 0 {k G Ng) are Hurwitz. This 
means that there exists a positive real number ei such that Re[crfc^i] < —ei(< 0) 
for all fc G Nq. In the limit of fc — )■ oo, the eigenvalue Afc —)■ oo, and n — 1 roots of 
PXk{s) converge to those of n{s) = 0 and one root converges to — oo, known as 
Butterworth pattern 18 , since n{s) is the n — 1-th order polynomial as shown in 
Lemma 1. Note that this rules out the possibility that the root asymptotically 
approaches to the imaginary axis. Let the roots of n(s) = 0 be {(Xao.if'iZi ■ 
Since n(s) = 0 does not have roots on the imaginary axis from Assumption 1, 
there exists a positive real nnmber e 2 such that Re[croo,i] < —e 2 (< 0) for all 
i = 1,2, • • • , n — 1. Substituting (16) into 0 and taking the integral over the 
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space leads to 


/ \\x{um= 

In Jn 


oo / n \ 

( Y Ck,ie‘"'‘-'*Wk{0) j (/)fc(0 





^ Wfc(0)(^fc(O 


fe =0 




i(^,0)|l d^, 


(17) 


where C is some constant and e := min(ei,e 2 )- Note that (j)k{0 = cos{kn^/Li) 


and the summation over i in (11) is dropped due to one dimensional spatial 
domain. 

(<^=) Suppose one of the polynomials p\^. (s) = 0 is not Hurwitz, which implies 
the existence of fcp S Nq and Zq € {1, 2, • • • , n} such that > 0. This im¬ 

plies that the growth rate o f the mode (pkg ) is non-negative, and the right-hand 
side of the first equality of (17) is bounded from below by Ce^®['^'=o-‘oh||u;^^(o)|| 
with some constant C. □ 


Proposition 1 allows us to analyze the stability of the infinite dimensional 
system ^ using the finite dimensional subsystems T,k/ {k = 0,1,2,-- - = 

0,1, 2, • It should be noted that there are still inhnite subsystems (fc = 
0,1, 2, • • • , ^ = 0,1, 2, • • •), although each subsystem is a finite dimensional sys¬ 
tem. In what follows, we show that the root locus method 18 allows us to 


systematically analyze the stability of the infinite number of subsystems 'Sk^e- 
It follows that the roots of p\^ ^ (s) = 0 (fc = 0,1,2, • • • , ^ = 0,1,2, • • •) lie on 
the root locus of the polynomial pa(s) = d(s)-l-An(s), which is the characteristic 
polynomial of a feedback system composed of h{s) and a feedback gain A (see 
Fig. i (Right)). Moreover, the definition implies that the feedback gain 
A = Xk^e monotonically increases in terms of k and i. Thus, the root locus 
starts from the point where A = Aq^o = 0 converges to the point where 
A = limfe_^_>oo Xk,e- This observation leads to the following graphical algorithm 
for stability analysis. 


Algorithm 1. Given the dynamics of a single cell h{s), we draw the root locus 
TZ of a negative feedback system composed ofh{s) and a constant feedback gain A 
with varying A S [0, oo]. If the root locus TZ lies inside the open left-half complex 
plane C_, the molecular communication network 0 is stable. 

This algorithm allows us to check the local stability of the molecular communi¬ 
cation network system ^ using the transfer function of a single cell h{s) and 
its root locus. Thus, the stability analysis of the infinite dimensional system is 
greatly simplified. 


Remark 2. The reason that the root locus gives only a sufficient con- 
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dition for stability is because the root locus is continuous in terms of the 
feedback gain A, while the poles of the subsystems T,k,e {k = 0 ,= 
0,1, 2, • • •) are discrete. In practice, however, the length of the spatial domain 
Li and L 2 is sufficiently large such that the discrete feedback gains = 

(fc = 0,1, 2, • • • , £ = 0,1, 2, • • •) are close to each other 
as shown in Section [3.2[ and the above algorithm can be practically considered 
as necessary and sufficient. 


It is worth noting that the root locus can be systematically drawn based on 
a set of known rules 18 . Specifically, the locus TZ starts from the n poles of 


h(s), and n — 1 roots converge to the zeros of h{s) and one pole goes to +00 or 
— 00 , since the relative degree of h{s) is one. 


Lemma 2. The root locus TZ starts from spec (A), and n—1 roots converge to 
spec(Acirc)- 


This is a direct consequence of Lemma 1 and the properties of the root locus 


18 


This lemma implies that the start and the terminal point of the root 
locus TZ can be characterized by the dynamics of the entire circuit and the 
internal circuit, respectively. This intuition is helpful from a synthetic biology 
viewpoint to determine the structure of the internal biocircuit to guarantee the 
stability/instability of the system. 


Remark 3. Arcak 22 showed algebraic conditions to guarantee the conver¬ 


gence of reaction-diffusion systems to spatially uniform solutions. When applied 
to our problem. Theorem 1 of 22 provides a sufficient condition for all the sub¬ 


systems with fc -I- £ > 1 to be exponentially stable. Thus, the condition is 
sufficient for the stability of the spatially homogeneous equilibrium x, provided 
that A is Hurwitz, or equivalently is stable. However, this sufficient condi¬ 
tion is conservative compared to Proposition 1, which is necessary and sufficient, 
and Algorithm 1 of this paper, since it is based on constituting a common Lya¬ 
punov solution that verifies the stability of all of the subsystems On the 

other hand, theorems shown in 22 can be used to certify the spatially unifor¬ 


mity of non-stationary solutions such as oscillations, while Proposition 1 cannot. 


3.2 Example: activator-repressor-diffuser system 

We illustrate the graphical stability analysis method using a hypothetical bio¬ 
chemical network shown in Fig. The system is composed of activator- 
repressor internal circuit and a single diffuser molecule that allows cell-to-cell 
communication. In Fig. A, ii and D stand for activator, repressor and dif¬ 
fuser, respectively. In the context of synthetic biocircuit design, one can choose 
an AraC protein as an (inducible) activator, i.e., A in Fig. and the araBAD 
operon as the corresponding promoter P^. For the repressor R, commonly used 
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> activation 
H repression 


I-1 t X 



I_ I T 



Figure 3: Schematic diagram of activator-repressor-diffuser system 


TetR and Lad proteins and the corresponding promoters can be used. The dif¬ 
fuser is implemented with AHL, which is a bacterial quorum sensing molecule. 


Let and denote the concentration of activator, repres¬ 

sor and diffuser, respectively. The dynamics of the molecular communication 
network system are then modeled by 


d = -5aa -I- 7a 
f = -SrT + Jr 


a?' Kr 

Kl + a2 A:2 + ^2 

Kl 

Kl + a^Kj + d^ 


d 


-Sdd + jd 


+ r 2 




(18) 

(19) 

( 20 ) 


where <5* {* = a, r,d) and 7 * (* = a, r,d) denote the degradation rates and the 
production rates of each molecular species, respectively. 


Case A. (Stable case): Suppose the half-lives of activator, repressor and 
diffuser are 19.8 minutes, 55.5 minutes and 138.6 minutes, respectively, which 
equates to 6a = 0.035min“^, 5r = 0.0125min“^ and 5d — 0.0050min“^. The 
production rates of each molecular species are ja = 2.5/iM • min“^, jr = 
1.65/J,M • min“^, jd = 0.25/j,M • min“^, and the Michaelis-Menten constants 
are Ka = Kr = Kd = lO/rM. The diffusion rate is fj, = 3.0 x 10“"^ mm^ • min“^. 
These parameter values are consistent with widely used parameters for numer¬ 
ical simulations in synthetic biology up to the order of magnitude (see for 
example). We leave the discussion on the realizability of such biomolecular com¬ 
munication networks in Remark 4 and focus on the numerical illustration of the 
root locus based analysis method here. 


Substituting the parameters into the model (18)-(20), we first calculate 
the spatially homogeneous equilibrium point of the system as [a*,r*,d*]^ = 
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[9.33,16.0,16.8]^. The Jacobian matrix A around this equilibrium point is 


0.243 

-2.93 ; 

0 

2.29 

-1.25 1 

-1.76 

0 

-0.757 1 

-0.500 


(21) 


The transfer function of a single cell h(s) is then calculated from (15) as 


h(s) 


+ 1.01xl0-2s +6.43x10-4 
s3 + 1.51 X 10-2s2 + 5.60x 10-4s + 3.54x 10-6 ’ 


Figure]^ (Right) illustrates the root locus TZ drawn by Algorithm 1. The 
eigenvalues of A and Adrc are spec(A) = {—0.00402 ± 0.0221j, —0.00702} and 
spec(Acirc) = {—0.00503 ± 0.0248}}, and we can see that these eigenvalues 
correspond to the initial and the terminal points of the root locus TZ. Since 
the root locus TZ stays in the left-half complex plane, the equilibrium point of 
the molecular communication network is locally stable. 

Figure]^ (Left) is the concentration gradient of the activator a{^,t) at t = 
2860 min. As expected from the graphical stability analysis result, the concen¬ 
tration converged to the spatially homogeneous equilibrium point. In the simu¬ 
lation, we used Li = 3mm and L 2 = 2mm, and the initial values a(4, 0), r(^, 0) 
and d(^,0) were set as a:*(l -I- 0.01 '0fc_^(^i, ^ 2 )), where x* is the 

equilibrium point of a(^, 0 ), r( 4 , 0 ) and (^, 0 ), respectively, and 


A Ail >^2)--cos y ) 



The full temporal dynamics is available in a movie format 33 , where the interval 
between frames is 5 minutes. We used COMSOL@ Multiphysics 4.4 for the 
simulation. 


Case B. (Unstable case): Let 7 ^ = 0.30/iM-min-4 and the rest of the param¬ 
eters being kept the same as Case A. The equilibrium point is then [a*, r*, = 

[7.63,15.6,14.5]^. The eigenvalues of A and Adrc are calculated as spec(A) = 
{-0.000169 ± 0.0237},-0.00792} and spec(Acirc) = {-0.00163 ± 0.0258}}. As 
shown in Fig. (Right), the root locus lies in the open right-half complex plane. 
In fact, multiple roots of p\{s) = 0 lie in the right-half complex plane, implying 
that the corresponding subsystems, or spatial modes, are unstable. 


Figure]^ (Left) illustrates the concentration gradient of the activator A at 
t = 2860 min. Since the equilibrium point is not stable, the concentration does 
not converge to the spatially homogeneous equilibrium. In fact, a coordinated 
time-varying spatial pattern appeared in this example, which implies that the 
system fell into limit cycle instead of converging to a spatially inhomogeneous 
equilibrium point. The full temporal dynamics is available in a movie format 
33 , where the interval between frames is 5 minutes. We used COMSOL(r) 


Multiphysics 4.4 for the simulation. 
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Figure 4: Case A: (Left) Concentration gradient of the activator A over cell 
population at t = 2860 min. (Right) Root locus 71 of the single cell dynamics 
h{s). 



Figure 5: Case B: (Left) Concentration gradient of the activator A over cell 
population at t = 2860 min. (Right) Root locus TZ of the single cell dynamics 
h{s). 


As Turing 15 pointed out, this coordinated pattern formation can be ex¬ 


plained by the fact that the subsystems corresponding to the non-zero spatial 
modes are destabilized. In the next section, we provide more rigorous argument 
of why and when the pattern formation can be expected. 


Remark 4. The libraries of activators 
ing sites 


34 


repressors 


35 and ribosome bind- 


36 recently became available for the use in synthetic biology. These 


libraries help us tune the rate parameters in a relatively predicted manner and 
would be useful for the implementation of the activator-repressor-diffuser bio¬ 
circuit. In fact, it is encouraging that multiple biocircuits were implemented 
and tuned using these parts libraries 34 35 37 . In the context of biomolecular 


communication network, the use of zinc finger proteins and small RNAs were 


also proposed in 38 . 


4 Conditions for Turing Pattern Formation 

4.1 Finite mode Turing instability 

As seen in the previous example, molecular communication networks can exhibit 
coordinated spatial patterns based on the passive communication mechanism. 
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Figure 6 : (Left) Concentration gradient of the activator A. (Right) Root locus 
TZ of the cell dynamics h{s). 


or diffusion. This is because non-zero finite spatial modes are destabilized and 
the growth rate of in ( 11 ) is positive. 

The instability of the non-zero spatial mode caused by diffusion is called Tur¬ 
ing instability named after Alan Turing, who first pointed out the mathematical 
basis of the emergence of inhomogeneous spatial patterns in reaction-diffusion 
In what follows, we introduce a formal definition of finite mode 


15 


equations 

Turing instability for biomolecular communication networks, then we investigate 
conditions for the Turing instability. 


Definition 1. (Finite mode Turing instability) Suppose Assumption 1 
holds. A biomolecular communication network 0 is finite mode Turing unstable 
around an equilibrium x, if 
0) PAo,o(s) = Po{s) = d{s) ^ 0; Vs e Co+, and 

(ii) Pa (7 -f s) = 0; 3A G (0, 00 ) and 3s G C+, where 7 = max(Re[so], 0) and sq 
is a root of n{s) = 0 with the largest real part. 

The condition (i) implies that the dynamics of each cell h{s) is stable around 
an equilibrium point x when there is no diffusion. The condition (ii) guarantees 
that the root locus of h{s) goes into the open right-half plane of the complex 
plane for some feedback gain A > 0 (see Fig. (Right)). In particular, the 
roots of n(s) = 0 are the terminal points of the root locus (Lemma 2), thus the 
condition (ii) also implies that the right-most roots of p\{s) = 0 are given when 
A is a non-zero and finite value. This guarantees that the dominant unstable 
spatial mode is finite and the resulting spatial pattern is of a non-zero finite 
spatial mode. We note that finite mode Turing instability is determined by 
the properties of the internal circuit Adrc and A but not the diffusion rate /i 
as shown in Definition 1. In practice, the wavelength of the dominant spatial 
mode needs to be longer than the size of a single cell. This can be tuned by 
varying the diffusion rate p after the finite Turing instability is guaranteed by 
the biocircuit structures Acirc and A. 

It is important to guarantee that the dominant spatial mode is finite, since it 
is physically and biologically implausible to achieve coordinated spatial patterns 
of infinite spatial mode at steady state. The following toy example shows the 
case where the right-most roots of pa(s) = 0 is given in the limit of A —)■ 00 . 
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Example: Suppose 



1.954 

-3.114 

; 0 

A = 

1.0 

-0.9540 

1 -1.223 


0 

1.0 

1 - 2.0 

3 anc 

L 2 = 2 . 

We can 

calculate 


(22) 


cell h{s) using A and the sub-matrix AcUc and obtain the root locus TZ of h{s). 
Figurej^ (Right) shows the right-most pole of the root locus TZ. As k,i ^ oo, or 
equivalently A —>■ oo, the root of px{s) = 0 converges to the root of n(s) = 0, or 
Sq. This implies that the growth rate of the state x is the largest for the infinite 
spatial mode, which is not plausible in real physical and biological systems. In 
fact, the simulation result in Fig. (Left) shows spatial patterns with non- 
realistic high spatial frequency. The time axis of the simulation, t, was scaled 
by r = 0.005t. The full temporal dynamics is available in a movie format 
We used COMSOL(r) Multiphysics 4.4 for the simulation. 


33 


It should be noted that the importance of the finite mode instability was 
not actively discussed in many classical works 14 39 where the systems are 


composed of n = 2 molecular species and both can diffuse in the spatial domain. 
This is because the dominant spatial mode is always guaranteed to be finite 
mode, no matter how the reaction rates are varied. 


4.2 Condition for finite mode Turing instability 

In this section, we show that at least three molecular species, Mi, M 2 and M 3 , 
are necessary to achieve finite mode Turing instability in molecular communica¬ 
tion networks. This implies that the well-known activator-repressor model [Tl] 
fails to make spatial patterns when there is only one diffusible molecule in the 
system. 

Assumption 2. The matrix A in the model is Hurwitz. 

This means that the reaction dynamics of a single cell h{s) is stable. As 
Definition 1 shows, the stability of h{s) is a necessary condition for finite mode 
Turing instability, and the existence of spatial patterns is essentially determined 
by the condition (ii) of Definition 1. Hence, by imposing Assumption 2, we 
hereafter focus on the condition (ii) of Definition 1. 

The following lemma provides the necessary and sufficient condition for finite 
mode Turing instability. 

Lemma 3. Consider the biomolecular communication network system and 
its linearized system Suppose Assumptions 1 and 2 hold. Then, the system 
0 is finite mode Turing unstable around an equilibrium x, if and only if either 
of (a) or (b) holds. 
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(a) Acirc is Hurwitz, and Re[pA(j<j-’)] = 0, Im[pA(jw)] = 0 for some A > 0 and 
some w S K. 

(b) ^circ is not Hurwitz, and Re[pA (7 + iw)] = 0, Iin[pA (7 + jw)] = 0 for some 
A > 0 and some a; S K, 7 ;= maXj,gspec(Aci„) 

This lemma is a direct consequence of Definition 1. The conditions (a) and (b) 
are equivalent to (ii) of Definition 1 with 7 = Re[so] and 7 = 0, respectively. 

Using this lemma, we can derive the following theorem showing that at least 
three molecular species are necessary to achieve hnite mode Turing instability. 

Theorem 1. Consider the biomolecular communication network system 0 
and its linearized system Suppose Assumptions 1 and 2 hold. Then, the 
system 0 composed of n = 1 and n = 2 molecular species cannot he finite mode 
Turing unstable. 


The proof of this theorem can be found in Appendix. This theorem implies 
that the activator-repressor-diffuser system illustrated in Section 3.2 is one of 
the minimal biomolecular communication networks that can exhibit finite mode 
Turing instability with n = 3. 


Moreover, the following theorem shows that hnite Turing instability is always 
caused by Hopf-Turing bifurcation, implying that not only coordinated spatial 
patterns but also temporal oscillations of concentration gradients are expected 
at steady state. 


Theorem 2. Consider the biomolecular communication network system 0 
and its linearized system 0 - Suppose the system 0 is finite mode Turing 
unstable around an equilibrium x. Then, the right-most pole of the characteristic 
polynomial p\{s) = 0 is a pair of complex conjugate. 


Proof. We prove by contradiction. Suppose the right-most root of p\{s) = 0 is 
a real number cr, and it is given when A = Aq. That is, px^lpj) = 0. It follows 
from (i) of Dehnition 1 that all of the poles of h{s) he in C_. Moreover, (ii) of 
Dehnition 1 implies that cr > max(Re[so], 0) > 0, where Sq is a root of n(s) = 0 
with the largest real part. Therefore, a is greater than any real poles and zeros 
of h{s). This, however, contradicts with the property of the root locus that 
the locus exists on real axis to the left of an odd number of poles and zeros of 
h{s) [^. □ 

In the next section, we turn our attention to the synthesis of the minimal 
molecular communication networks with n = 3 molecular species and investigate 
their properties in detail. 
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5 APPLICATION TO THREE-COMPONENT 
CIRCUIT SYNTHESIS 

In this section, we consider molecular communication networks consisting of 
two internal circuit molecules, A4i and AI21 and one diffuser AI3. For this 
class of systems, we first derive an analytic condition for finite mode Turing 
instability. Using the analytic condition, we characterize a structural property 
of biomolecular communication networks that admit Turing instability. 

Let ai and j3i denote the coefficients of the following characteristic polyno¬ 
mials. 

ls /3 — A| = -f Ot2S^ -\- OiiS -f Og, (23) 

ls /2 — ^circl = s'^ + PlS + Po- (24) 

Then, a necessary and sufficient condition for finite mode Turing instability is 
given by the following theorem. 

Theorem 3. Consider the biomolecular communication network system 0 
composed of n = 3 molecular species and its linearized system 0. Suppose 
Assumptions 1 and 2. The system 0 is finite mode Turing unstable around an 
equilibrium x, if and only if the following (A) or (B) is satisfied. 

(A) /?! > 0, /3o > 0, 

+ Pia 2 — /3o < —2\/ Pi{aia 2 — Og), 

(B) /?! < 0, - 4/3o < 0, 

~ Pi fdo A- PiOi2 — Qfi > 0. 


The conditions (A) and (B) in Theorem 3 are equivalent to (a) and (b) of 
Lemma 3, respectively. The complete proof is shown in Appendix. It should be 
noticed that Pi {i = 0,1) are determined from the parameters of the internal 
circuit Acirc) while ai (i = 0,1,2) depend on the entire system including the 
diffuser. Theorem 3 allows us to further investigate the properties of three- 
component systems. In what follows, we show that a molecular communication 
network needs to have a certain structural property to admit finite mode Turing 
instability. 

We consider a class of cascaded three-component circuits illustrated in Fig. 

The matrix A is then defined by 


«11 

012 

0 

021 

022 

O23 

0 

O32 

^diff 


(25) 


The activator-repressor-diffuser network in Section 3.2 is a class of the cas¬ 
caded circuit with 012 < 0, 021 > 0, 023 < 0 and 032 < 0. Note that the signs 
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Figure 7: Schematic diagram of three-component circuit 


of 012,021,023 and 032 depend on whether the molecule and 7W3 are 

activator or repressor. The following proposition specifies the combinations of 
activator and repressors that admit finite mode Turing instability. 

Proposition 2 . Consider the biomoleeular communication networks with the 
cascaded structure shown in Fig. Suppose Assumptions 1 and 2 hold. Let 
ri := 012O21 and r 2 '.= 023O32. Then, the network admits finite mode Turing 
instability only if ri <0 and r2 > 0. 


The proof can be found in Appendix. This proposition narrows down the 
design space of biomoleeular communication networks. In other words, the 
combination of the molecules Mi, M 2 and At3 needs to be either of (acti¬ 
vator, repressor, repressor) or (repressor, activator, activator) in order for the 
system to achieve finite mode Turing instability. It should be noticed that the 
example shown in Section [3.2| is an example of (activator, repressor, repressor) 
configuration. 

6 Conclusion 

In this paper, we have proposed a control theoretic modeling and analysis frame¬ 
work for biomoleeular communication networks modeled by reaction-diffusion 
equations with a single diffusible molecule. We have first shown that the anal¬ 
ysis of the infinite dimensional system can be simplified by decomposing it 
into finite dimensional subsystems. We have then provided a graphical stabil¬ 
ity analysis tool using the root locus approach and demonstrated on the novel 
activator-repressor-diffuser system, which allows for the formation of inhomo¬ 
geneous concentration gradient over the cell population. The latter half of the 
paper has been devoted to analyzing the conditions for finite mode Turing in¬ 
stability. Using the analytic condition of Turing instability, we have shown that 
the network needs to have at least three molecular species to admit Turing in¬ 
stability. Moreover, the in-depth study of three-component circuits provided 
a necessary activation-repression network structure condition for Turing insta¬ 
bility. Our future work will be devoted to addressing robustness issues such 
as cell-to-cell variability and parameter uncertainties of biomoleeular communi¬ 
cation networks. The activator-repressor-diffuser system will also need further 
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study to characterize the parameter regions for Turing instability. 
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.1 Proof of Lemma 1 

Using the definition Q, we can write h{s) as 


h{s) = 


'’adj(s/„ - A)e„ 
\sln - 71 ] 


(26) 


where adj(A) denotes the adjugate matrix of X. It follows from the definition 
of the adjugate matrix that the (n, n)-th entry of adj(s/„ — A) is — Acird, 

and thus the numerator of h{s) can be written as e^adj(s/„ — A)e„ = — 

Acircj- n 


.2 Proof of Theorem 1 

For n = 1, it follows that the molecular communication network is finite mode 
Turing unstable if and only if the characteristic polynomial px{s) = s — Odiff + A 
satisfies (a) of Lemma 3. We can easily see that Re[pA(/w)] = 0 cannot be 
satisfied for all A > 0 and a; G K, since Odiff < 0. 

For n = 2, we define coefficient Oq, ai and do by ls/2 — A] = -I- Ois -I- oq 

and ls/2 — Acircj = s + do- In what follows, we show that both (a) and (b) of 
Lemma 3 cannot be satisfied. 

Case 1: Adrc is Hurwitz: It follows that Im[pA(/w)] = uj{ai -I- A) for all 
A > 0 and w ^ 0, since Oi > 0 holds when A is Hurwtiz. For w = 0, Re[pA(0)] = 
ao -I- Ado 0 holds for all A > 0 because of oq > 0 and do > 0. Thus, the 

condition (a) in Lemma 3 cannot be satisfied. 

Case 2: Adrc is not Hurwitz: It follows that Im[pA(7+jA>)] = uj{2^ + ai+\), 
where 7 = do(> 0) from the definition. Thus, Im[(A,7-I-jw)] 0 for all 

A > 0 and w 0, since oi > 0 and 7 > 0. When w = 0, Re[pA(7)] = 

7^ -I- 7(01 -I- A) -I- oq -b Ado > 0 hold, and Re[pA(7)] y^ 0 hold for all A > 0. This 
completes the proof. □ 


.3 Proof of Theorem 3 

We prove the theorem by calculating conditions (a) and (b) of Lemma 3 for 
n = 3. 
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In what follows, we derive the conditions (A) and (B) by considering the 
case where Acirc is Hurwitz and not Hurwitz, respectively. We can easily verify 
that neither (a) nor (b) of Lemma 3 can be satisfied when w = 0. Hence, we 
hereafter show the proof for the case of w 7 ^: 0 . 

Case 1: Adrc is Hurwitz: A necessary condition for Adrc being Hurwitz is 


/3i > 0 and /?o > 0 . 

Moreover, px{joj) = 0 implies that 

Re[px{ju;)] = -a 2 W^ + ao + A(-a;^ + /Jg) = 0 
lm.[pxijuj)] = -a;{a;^ - (ai + A/3i)} = 0. 


(27) 


(28) 

(29) 


It follows from (29) that = ai+A/3i(> 0), since w 7 ^ 0. We can then eliminate 
from (28), and we have the following quadratic equation of A. 


/3iA^ + (ai + (3ice2 — /3o)A + q:iq ;2 — 0:0 — 0. 


(30) 


Since /3i > 0 holds from (1^, and aia 2 — ao > 0 holds from Assumption 2, 


a necessary and sufficient condition for (30) having a positive real solution is 
given by the following (Cl) and (C2). 

(Cl) The determinant of (30) is non-negative. That is, 

((ai -f Pi)a2 - PlY - AI5i{aia2 - ao) > 0. 


(31) 


(C2) The coefficient of A in (30) is negative. That is, 

ai + Pia2 - /3o < 0. 


Summarizing (31) and (32), we have 


ai + (3ia2 — /3o < —2\/ /3i{aia2 — ag). 


(32) 


(33) 


The condition (A) is obtained from (27) and (33). 


Case 2: Adrc is not Hurwitz A necessary condition for Adrc not being 
Hurwitz is 


/3i < 0 or /3o < 0. (34) 

Let 7 ± jp ( 7 , p > 0) denote a pair of eigenvalues of Adrc with the largest real 
part. Then, [(7 + jp)l 2 — Acird = 0 implies 

7^ + di7 + Po-p‘^ =0 (35) 

(27 + /3i)7 = 0. (36) 
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In what follows, we first show that pa(7 + jv) 0 for all A > 0 and uj when 
r] = 0. It follows from Im[pA(7 + jw)] = 0 that 

= 37 ^ + 2a27 + Q!i + A(27 +/3i). (37) 

Substituting this into Re[pA(7 + j^)] = 0 yields the following equation of A. 

(27 + /3i)A^ + { 97 ^ + (402 + 3/3i) 7 + oi + /3ia2}A 


+ {87^ + 8027^ + 2(ai + 02)7 + (aia2 - ao)} = 0. 


( 38 ) 


It should be noted that the coefficient of A^ is non-negative, or 2j + Pi > 0, since 
7 is the largest real root of -|- /3is -I- /3o = 0 from the definition. Moreover, 
the coefficients of A and the constant terms of (38) are also positive because 
oo > 0, 02 > 0 and 0102 — oq > 0 from Assumption 2, and (34) holds. Thus, 
Pa (7 + jv) 7 ^ 0 for all A > 0 and w when r] = 0. 

When rj ^ 0, the determinant of |s/2 — Adrcl = 0 is negative, thus we have 

Pf - 4/3o < 0. (39) 

Moreover, Im[pA(7 + jw)] = 0 implies 

07^=37^-1-202/3-1-01. (40) 

Substituting this into Re[pA(7 + jw)] = 0 yields 


{ 7 ^-( 3 , 8 ^ -I- 202/3 -I- oi)}A - {87^ -I- 8027^-1- 
2 (ai -I- 02)7 -I- (ai02 — oo)} = 0 . 


( 41 ) 


Note that the constant term of ( 41 ) is negative because of Assumption 2 and 
7 > 0 . Thus, ( 41 ) has a positive real root if and only if 

rf' — (37^ -I- 2o 27 + ai) > 0 - ( 42 ) 

□ 


The condition (B) is obtained from ( 34 ), ( 35 ), ( 36 ), ( 39 ) and ( 42 ). 


.4 Proof of Proposition 2 

We first show ri < 0 . Suppose (A) of Theorem 3 holds. It follows from / 3 i = 
— (ail + 022) > 0 of (A) that either or both of an and a22 need to be negative. 
It follows that ri = 012021 < 0 when on and O22 have opposite signs, hence we 
consider the case of an < 0 and a22 < 0 in the following. 

We first note that 0102 — Oq > 0 holds from Assumption 2 , and it can be 
equivalently written as 

(on + 022)^1 -I- (022 + adiff)r2 

— (on -I- a22)(a22 + Odiff)(adiff + an) > 0 . ( 43 ) 
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Dividing (43) by an + 022 (< 0), we have 


ri < (022 + adiff)(a(iiff + an)-^-7’2. (44) 

an + a22 

On the other hand, r 2 satisfies the following inequality from ai + idia 2 — /3o < 
-2y//3i(aia2 - ao) < 0 of (A). 

(an + a22)(an + 0,22 + 2 adiff) - r2 < 0 . ( 45 ) 

Since an < 0,a22 < 0 and Odis < 0, the second term on the right side of ( [44| 
satisfies 


022 + Odiff 
On + 022 


1"2 < —(022 + adiff)(an + 022 + 2 adiff)- 


(46) 


Summerizing (44) and (46), we have 


ri < (022 + adiff)(odiff + an) - 

On + 022 

< (022 + Odiff)(adiff + an) 

— (022 + adiff)(an + 022 + 2 adiff) 

= ~(022 + Odiff)^ < 0. 


(47) 


Suppose (II-B) of Theorem 3 holds. It then follows from — 4/3o < 0 of (B) 


that 


(an + 022)^ - 4(011022 - n ) = (an - 022)^ + 4 ri < 0 . 

This concludes ri < 0. 

Next, we show that r 2 > 0 is necessary. Suppose (B) of Theorem 3 holds. 
Then, the direct calculation leads to —/3i + ,do + /3ia2 — ai = r 2 > 0. In what 
follows, we consider the case where (A) of Theorem 3 holds. 

It follows from adiff < 0 and /3i = —(ai+a 2 ) > 0 of (A) that (ai +a 9 + 2a.i) < 
0. Then, ai + /3ia2 — /3o < —2-\//3i(aia2 — ag) < 0 of (A) implies (45). Thus, 
^2 > (oi + a2)(ai + a2 + 2a3) >0. D 
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